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Wave propagation in the magnetic sun 

By T. Hartlep, M. S. Mieschf, and N. N. Mansour 



This paper reports on efforts to simulate wave propagation in the solar interior. Pre- 
sented is work on extending a numerical code for constant entropy acoustic waves in 
the absence of magnetic fields to the case where magnetic fields are present. A set of 
linearized magnetohydrodynamic (MHD) perturbation equations has been derived and 
implemented. 



1. Introduction 

The evolution of the solar interior is actively studied theoretically, numerically, and 
observationally. Many recent advancements in our understanding of the structure and 
dynamics of the sun have been made by studying solar oscillation. This is the science of 
helioseismology. Through ground- and satellite-based telescopes, oscillations on the solar 
surface are being observed and, through appropriate techniques, are used to infer infor- 
mation about the sun's internal structure a nd composition. One in teresting technique 
provides so-called far-side images of the sun ( Lindsev fc Braun 200Ci) . where maps of ac- 
tive regions on the far side (the side of the sun facing away from earth) are inferred from 
the oscillations observed on the front side (the one facing earth). For a review of other 
helioseismic methods, as well as general properties of the observed solar oscillations, see 
Christensen-Dalsgaard In general, these inferences are based on simplified mod 



els of wave propagation, and researchers in the field fe.g.. IWerne. Birch fc JulienI 120041 ) 
have expressed their need for wave propagation simulations to test and calibrate these 
methods. 

Of course, wave propagation can be studied by simulating the full compressible equa- 
tio ns; for instance, simulatio ns of the shallow upper layer of the solar convection zone 
bv IStein fc Nordlund ( 2000l ) demonstrate excellent agreement with existing analytical 
theories and observations. But due to the enormous computational requirements, these 
types of simulations are able to simulate only a very small part of the sun. For global sim- 
ulations such an approach is not expected to be feasible in the immediate future. Except 
for the very near surface region, though, flow velocities in the sun are much smaller than 
the speed of sound and a perturbation approach is applicable. Thus, acoustic waves can 
be treated here as small perturbations traveling through a base flow. This is the foun- 
dational idea of the this project. The base flow itself can be either artificially prescribed 
or it can be provided by other simulations such as three-dimensional global simulations 
of solar convection in the ane lastic approximation ( Miesch 1998 : Brun fc Toomrd 2002t 
Brun. Miesch fc Toonir3l2004[ ). Our work extends that of iHartlep fc Mansouil (|2005l i^ 
including the effects of a prescribed magnetic field on the propagation of waves, which 
is especially of interest for testing the far-side imaging technique previously mentioned. 
We present here the perturbation equations we derived and our numerical method, and 
comment on the current status of this project. 
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2. Model equations 

The idea underlying the derivation is to spht the dependent variables (density, pressure, 
magnetic field strength, etc.) into base variables, denoted by ~, and wave perturbations, 
•', and then to derive appropriate equations for the perturbations. We will take the base 
state to be non-rotating, without flows (v = 0), and with prescribed three-dimensional 
sound speed c^, density p, magnetic field B, and electric current distribution J . The real 
sun, of course, is rotating and has non-vanishing flow velocities. We plan to add these 
influences at a later stage, but the current focus is on the effects that magnetic fields and 
spatial variations in the speed of sound have on the propagation of acoustic waves. The 
linearized perturbation equations arising from the continuity and Euler's equation are: 

V-m' (2.1) 
Vp' + p'g + L' (2.2) 

with m' = pv' , g = —rg, and L' being the mass flux, the acceleration due to gravity, 
and the perturbation of the Lorcntz force, respectively. The background state is assumed 
to be in magnetohydrostatic balance: Vp = —pg + L. In principle, another term in the 
momentum equation (|2.2[) appears to account for the perturbation of the gravitational 
acceleration caused by the waves, ^pg'. This term can be expected to be very small, 
though, and is neglected here. Needed next is an equation that relates the pressure to 
the density perturbation. Assuming adiabatic processes, we have p' = cj.p' + jp/Cps' 
with Cp and 7 being the specific heat at constant pressure and the adiabatic index, 
respectively. Only isentropic waves are considered at this point, and therefore s', the 
entropy perturbation, is set to zero. The magnetic field, of course, is divergence free, 

V ■ J5' = 0, (2.3) 



dt 



and is given by Faraday's law of induction: 

^B' = -cV X E'. (2.4) 

dt ^ ' 

Also included are Ohm's law (in the infinite conductivity limit), Ampere's law, and the 
Lorentz force given by: 

E' = X B, (2.5) 

c 

J' = X B', (2.6) 

L' ^-\J' X B + J X B']. (2.7) 
c ^ ^ 

The displacement current dtE' in the equation for the electric current has been ne- 
glected because the velocities involved are small compared to the speed of light c. These 
linearized equations do not account for the initial excitation of the waves, which is a 
non-linear process. Acoustic waves are expected to be generated by the vigorous turbu- 
lent convection close below the solar surface. We simulate these stochastic excitations by 
artificially adding a random forcing term to Eq. (|2.ip . 
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3. Numerical method 

The preceding equations are solved in spherical coordinates with a pseudo-spectral 
method. Scalar quantities such as pressure and density are expanded in t erms of spher- 
ical harmonic basis functions for their angular structure, and B-splines (Ide Boon Il987 ; 
Loulou. Moser. Mansour fc Cantwelilll997t iHartlep fc Mansouilliool l2005l) for their ra- 



dial dependence. Vector fields such as the magnetic field, the Lorentz force, and the mass 
flux are expanded in vector spherical harmonics and B-splines. Vector spherical harmon- 
ics were selected because of the coordinate singularities in spherical coordinates, which 
are most easily treated in those basis functions (see Appendix Typically, B-splines of 
polynomial order 4 are used, and the spacing of the generating knot points is chosen to 
be proportional to the local speed of sound. This results in higher radial resolution near 
the solar surface, where the sound speed is low (less than 7 km/s) compared to the deep 
interior, where the sound speed surpasses 500 km/s. The numerical domain is chosen to 
be larger than the solar radius of approximately 696 Mm by an additional 2-4 Mm to 
account for an absorption layer as a means of realizing non-reflecting boundary condi- 
tions. This layer is implemented by adding damping terms —ap' and —am' to Eq. (|2.ip 
and p.2p . where cr is a positive damping coefficient that is independent of time, which is 
set to be zero in the interior, and from the solar surface up increases smoothly into the 
buffer layer. 

The equations are then recast using an integrating factor exp{at), and advanced using 
a staggered leapfrog scheme in which p' and B' are advanced at the same time, and m' 
is offset by half of a time step. The basic properties of the background state, p, g, and c^, 
are obtained fro m solar models, which provide s hell-a. veraged properties of the quiet sun. 
Solar model S bv lChristensen-Dalsgaard et al. ( 1996f) is used for radii below 696.5 Mm, 
the highest radius consider ed in that mode l. Prop erties above this radius are taken from 



chromosphere model C by IVernazza et 



(|l98lh . Models for 3-D spatial variations of 
the temperature, and therefore the sound speed, and the magnetic field, in structures 
such as sunspots are then added to this basic, spherically symmetric background. 



4. Summary 

We have derived a set of linearized MHD perturbation equations for the propagation 
of waves in a simplified sun that is non-rotating and has no background fiows, and we 
have implemented these equations as an ext ension to our previous non-magnetic wave 
propagation code (jHartlep fc Mansouii 120051 ) . This represents a considerable extension 
and required developing routines for vector spherical harmonics, including fast trans- 
formations between spectral and physical space (which were not present in the original 
code). The new implementation is still in the testing phase, but we hope to be able to 
perform numerical experiments (e.g., a test of the far-side imaging technique) soon. 



Appendix A. Treatment of coordinate singularities in spherical coordinates 

The sun is very well approximated by a sphere, and spherical coordinates are therefore 
the obvious choice in the numerical method. This requires, however, that special care is 
taken to ensure regularity of physical quantities on the polar axis and at the coordinate 
center. The use of a spherical harmonics expansion guarantees smoothness at the polar 
axis, but additional constraints on the expansion arise from requiring regularity at the 
origin. Fortunately, in many applications the origin can be excluded from the numerical 
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domain and therefore those constraints are irrelevant. In the present case though, acoustic 
waves travel through the whole sphere, and the center cannot be excluded, requiring that 
we enforce these constraints. The constraints are different for scalars such as density and 
pressure, and for vector quantities such as velocity and magnetic field. The following 
briefly outlines their derivation. 

We start with a scalar variable, which is expanded in terms of spherical harmonics 
Yi „i for its angular dependence and, for each spherical degree I and azimuthal index m, 
some radial functions fi^m'- 

oc +/ 

Of course, in the actual implementation the infinite sum is truncated at a chosen maxi- 
mum value of 1. Using the definition of the spherical harmonics and the transformation 
to Cartesian coordinates [x = r sin cost/), y — rsin^sin^, z = rcosO) we can rewrite 
Eq. (jAll) into: 

OO +/ / 

^— m—~l k—kmin 

with some constants Ck,i,m, and where kmin = (/ + |to|)/2 for even values of / + \m\, and 
kmin = + |m| + l)/2 for odd values of Z + |m|. The plus and minus signs in Eq. (jA 2^ 
are used for positive and negative values of m, respectively. The terms involving x, y, 
and z are always regular at the center since |m| and 2fc — (/ + |m|) are non-negative. The 
only constraint is that r~'^^~^^ fi^rn{f) must also be regular. This function, expanded in a 
Taylor series around r = 0, reads as: 

oo 

r-^''+' krn{r) - <^i.m,pr-^''^'r^, (A 3) 

where individual terms are only smooth at the origin if the combined exponent —2k + l+p 
is positive and even. For all terms for which this is not the case, the expansion coefficient 
o^i,m,p must vanish. In other words, the constraint is that the radial function fi^m must 
be of the form: 

fiMr)=r'P{r'), (A 4) 

with P(r^) bein g a smooth fun ction in r^. An alternative derivation of this constraint 



can be found in iStanawavl (|198S). For similar c onstr aints that arise in Fourier expansions 



in cylindrical coordinates see lLewis fc Bellan 



Additional complications arise for vector quantities from the choice of unit vectors. A 
rather simplistic way would be to write a vector field in terms of r-, 9-, and i/i-components, 
and then to expand these components like scalars individually in spherical harmonic 
functions, i.e., 

oc +/ 

F{r,0,cP)=Yl E i1.™(^,</'){/r,z,m(r)r(0,0)-f/e,z,™(r)0(0,0)+/<^,i,™(r)0(0, ,/))}. 

1=0 m=-l 

(A5) 

Such an expansion is fine if the computational domain does no t include th e origi n r = 0, 
as is the case in the simulations of the solar convection layer bv lBrun et al. ( 2004 ). In our 



case, however, it is not an appropriate choice, since, for the vector field to be regular at 
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the origin, fr,i,m, fed,m, and f^^i^m are not independent of each other and the regularity 
condition becomes hard to enforce. The alternative here is to use so-called vector spherical 
harmonics as angular basis function: 

F(r,0,0) =^ {^xxn^{r)Xu„,{e,c^)+fvxn^{r)Vu^{e,cj,)+fwxn^{r)Wu^{e,<p)], 

1=0 m=-l 

(A 6) 

where Xi^m, and Wi^m are defined according to iHilll (|l954l ) by 



1 mYi,r^{0,cf,) ^ 1 .d^^ 



+ sine + 80 



7 / 1 tmYi^rajO,^) 

{I + 1){21 + 1) sine ' ^ 



+ ^V ^(2^ + 1) sing ■ 

Here, the regularity conditions for the radial functions decouple and are actually very 
similar to what is found for scalar fields. For instance, X; m can be rewritten as 

Xl^r,^{e,<t)) = [x + iy)]^y/{l + ni){l - m + l)Yi^m^i{e, 4>) 
1 



+ {x- iy) - y/{l-m)il + m + l)Yi^,n+i {9, 4>) 
+ iTOll,„,(0,0), (A 10) 

and therefore, following the arguments above, the radial function fx,i,m must behave 
just like a scalar, i.e., 

fx,Lm{r)^r'Px,i,m{r''), (All) 

with again Px,i,m{i'^) being a smooth function in r^. The results for the other two 
functions are: 

/vj.mW = r'+iPv,i,m(r2), (A 12) 

(r) = r'-^PwdAr^)- (A 13) 

In our numerical method, all radial functions are expanded in B-splines, which are piece- 
wise polynomials. Conditions (|A 4p and (jA lip - (|A 13|) result in a linear coupling of the 
expansion coefficients for B-splines near the center. These linear systems are very small 
and easy to implement. Only the first fc -f 1 coefficients are coupled, where k is the poly- 
nomial order of the B-splines used (here usually fc = 4) , since these are the only B-splines 
that are non-zero close to boundary r = 0. 
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